EMG - Overview
Difficulty Level:
Tags other☁overview☁emg

The action of skeletal muscles for movement and postural control is voluntary, causing drastic differences in EMG signal when comparing to ECG, namely the inexistence of natural periodicity.

For contracting a muscle , a large set of motor units needs to be activated, so that the acquired EMG signal is the sum of their elementary potential changes. Because of this "summation" process, EMG seems to be a little "anarchic", and the essence of EMG signal processing is in study the activation zones .

This Jupyter Notebook provides an overview on the processing steps to be taken into account for EMG signal processing for the extraction of parameters as well as event detection. It implies simple techniques, which can be extended with more complex processing methods.

The following topics will be covered:

A. Pre-Acquisition
B. Acquisition
C. Signal Loading
D. Pre-Processing
E. Event Detection
F. Parameter Extraction

Importing relevant packages

In [4]: 
# Biosignalsnotebooks python package
import biosignalsnotebooks as bsnb

# Scientific packages
from numpy import max, min, average, std, sum, sqrt, where, argmax, absolute, array, random, zeros
from scipy.integrate import cumtrapz
from scipy.signal import welch

A. Pre-Acquisition

A.1 - Placement of EMG electrodes for monitoring muscular activity: Adductor pollicis
The 1st ossa metacarpalia is the "reference" structure. The electrodes should be placed at the middle of the 1/4 distal region .

Adductor pollicis

A.2 - Selection of a bone structure away from the active muscles: Clavicle

Ground electrode on clavicle

B. Acquisition

B.1 - Sampling Rate:
According to the Nyquist Theorem , to ensure that each signal (EMG, ECG, EDA...) is acquired correctly (avoiding aliasing) the sampling rate should be at least the double of the maximum frequency component present in the signal, being this threshold known as "Nyquist Rate" .

In the following images, we can see the effect of sampling rate choice. EMG muscular activation was collected at 10 Hz and at 1000 Hz .

In [5]: 
# Hidden packages import.
from numpy import linspace

# [1000 Hz - 16 bit]
# Load of data
data_1000_hz = bsnb.load_signal("emg_1000_hz_16_bits")

# Identification of the mac-address of the device used to collect data through a 1000 Hz sampling rate.
mac_1000 = list(data_1000_hz.keys())[0]

# Storage of the label of the desired channel.
channel_1000 = list(data_1000_hz[mac_1000].keys())[0]

# Definition of the starting and end points of the visualisation window.
sample_rate_1000_hz = 1000
start_sample_1000_hz = int(sample_rate_1000_hz * 3.5) # 3.5 seconds
end_sample_1000_hz = int(sample_rate_1000_hz * 3.7) # 3.7 seconds

# Select the intended channel/window.
data_1000_hz = data_1000_hz[mac_1000][channel_1000][start_sample_1000_hz:end_sample_1000_hz]

# Time-axis.
time_1000_hz = linspace(0, len(data_1000_hz) / sample_rate_1000_hz, len(data_1000_hz))

# [10 Hz - 16 bit]
# Load of data
data_10_hz = bsnb.load_signal("emg_10_hz_16_bits")

# Storage of the label of the desired channel.
channel_10 = list(data_10_hz.keys())[0]

# Definition of the starting and end points of the visualisation window.
sample_rate_10_hz = 10
start_sample_10_hz = int(sample_rate_10_hz * 3.5) # 3.5 seconds
end_sample_10_hz = int(sample_rate_10_hz * 3.7) # 3.7 seconds

# Select the intended channel/window.
data_10_hz = data_10_hz[channel_10][start_sample_10_hz:end_sample_10_hz]

# Time-axis.
time_10_hz = linspace(0, len(data_10_hz) / sample_rate_10_hz, len(data_10_hz))

# Plotting data.
bsnb.plot([time_10_hz, time_1000_hz], [data_10_hz, data_1000_hz], grid_plot=True, grid_lines=1, grid_columns=2, y_axis_label="RAW Data", title=["Sampling Rate: 10 Hz", "Sampling Rate: 1000 Hz"])

Observation of different sampling frequencies:

  • 10 Hz : Aliasing due to low sampling rate - in the 0.2 s window exists a clear aliasing problem, taking into consideration that only 2 samples were acquired, which is not enough to properly represent the EMG signal
  • 1000 Hz : The random nature of the motor unit activation (responsible for the muscle contraction) is properly observed - enough samples collected

For more information on sampling rate see Jupyter Notebook entitled "Problems of low sampling rate (aliasing)"

B.2 - Resolution:

Resolution is another parameter that must be configured prior to acquisition. For PLUX acquisition systems, that support acquisitions with resolutions between 6 ( bitalino ) and 16 bits, a smaller resolution may not cause easily observable changes in EMG or ECG, as we will see in the following images by comparing the same signal when acquired with distinct resolutions.

EMG muscular activation was collected at 8 Bit and 16 Bit with a sampling frequency of 1000 Hz .

In [6]: 
# Load of data
data, header = bsnb.load_signal("emg_1000_hz_8_16_bits", get_header=True)


# Identification of mac addresses and channels
# Acquisition 8 bits (channel 5)
mac_address_8_bits = list(data.keys())[1]
channel_8_bits = list(data[mac_address_8_bits].keys())[0]
# Acquisition 16 bits (channel 6)
mac_address_16_bits = list(data.keys())[0]
channel_16_bits = list(data[mac_address_16_bits].keys())[0]

# Sampling rate
sr = header[mac_address_8_bits]["sampling rate"]
device = header[mac_address_8_bits]["device"]

# Transposition of data to an independent variable
# [8 bits]
data_8_bits = data[mac_address_8_bits][channel_8_bits]
# [16 bits]
data_16_bits = data[mac_address_16_bits][channel_16_bits]

#Unit Conversion
# [8 bits]
data_8_bits = bsnb.raw_to_phy(sensor="EMG", device=device, raw_signal=data_8_bits, resolution=8, option="mV")
# [16 bits]
data_16_bits = bsnb.raw_to_phy(sensor="EMG", device=device, raw_signal=data_16_bits, resolution=16, option="mV")


#Generate Time Axis
time = linspace(0, len(data_8_bits) / sr, len(data_8_bits))


#Plot - adjust the plot axis and zoom like above
bsnb.plot_compare_resolutions(time, data_8_bits, data_16_bits, sr, y_axis_label="Electric Tension (mV)",
                             time_start=7.75, time_end=8.0, y_range=[-0.6, 1.3])

Observation about different resolutions at 1000 Hz:
The signal shows a slight difference when comparing the same signal with 8 Bit , where the signal is slightly more angular and 16 Bit , where the signal is slightly more round.
However in some cases, such as temperature monitoring see Jupyter Notebook entitled "Resolution - The difference between smooth and abrupt variations" the signal is considerably affected when different resolutions are chosen.

C. Signal Loading

C.1 - Load data from signal samples library

In [7]: 
# Load of data
relative_file_path = "../../signal_samples/emg_1000_hz_16_bits_solo.h5"
data, header = bsnb.load(relative_file_path, get_header=True)
In [8]: 
from sty import fg, rs
print(fg(98,195,238) + "\033[1mHeader\033[0m" + fg.rs)
print(header)
print(fg(98,195,238) + "\033[1mData Structure\033[0m" + fg.rs)
print(data)

Header
{'channels': array([1]), 'comments': '', 'date': '2019-10-14', 'device': 'channeller', 'device connection': b'BTH00:07:80:D8:A8:82', 'device name': b'00:07:80:D8:A8:82', 'digital IO': array([0, 1]), 'firmware version': 775, 'resolution': array([16]), 'sampling rate': 1000, 'sync interval': 2, 'time': '16:52:39.448', 'sensor': [b'EMG'], 'column labels': {1: 'channel_1'}}
Data Structure
{'CH1': array([32960, 33014, 32967, ..., 32793, 32743, 32723], dtype=uint32)}

C.2 - Information from file Header

In [9]: 
#get information which is stored inside variables
ch = "CH1" # Channel
sr = header["sampling rate"] # Sampling rate
resolution = header["resolution"][0] # Resolution (number of available bits)
device = header["device"]

C.3 - Store the desired data in an individual variable

In [10]: 
#RAW DATA
signal = data[ch]

C.4 - Convert the RAW data to values with a physical meaning (in this case electric voltage | mV)

In [11]: 
# Signal Samples Conversion
signal_uv = bsnb.raw_to_phy("EMG", device, signal, resolution, option="mV") # Conversion to mV

C.5 - Generate a time-axis

In [12]: 
time = bsnb.generate_time(signal)
In [13]: 
bsnb.plot([time[::10]], [signal_uv[::10]], y_axis_label="Electric Tension (mV)")

D. Pre-Processing

The effect of a bandpass-filter on the signal (power spectrum) is shown in the following images.

D.1 - Presentation of the original signal power spectrum
With this approach it will be possible to understand the frequency content of the collected data and understand where it is located the informational band.

In [14]: 
# [Baseline Removal]
pre_pro_signal = signal - average(signal)

# Power spectrum
freq_axis_orig, power_spect_orig = bsnb.plotfft(pre_pro_signal, sr, doplot=False)
bsnb.plot(freq_axis_orig, power_spect_orig, x_axis_label="Frequency (Hz)", y_axis_label="Relative Intensity (a.u.)")

D.2 - Inclusion of some artificial Gaussian noise to the signal

In [15]: 
art_signal = signal + random.random(len(signal)) * 0.035 * max(signal)
In [16]: 
# [Baseline Removal]
art_signal = art_signal - average(art_signal)

# Power spectrum
freq_axis_art, power_spect_art = bsnb.plotfft(art_signal, sr)

# Generation of the two power spectra.
bsnb.plot([freq_axis_orig[::10], freq_axis_art[::10]], [power_spect_orig[::10], power_spect_art[::10]], x_axis_label="Frequency (Hz)", y_axis_label="Relative Intensity (a.u.)", legend_label=["Original Signal", "Signal Contaminated with Artificial Noise"], grid_plot=True, grid_lines=2, grid_columns=1)

D.3 - Application of a band-pass filter to remove high-frequency noise-related components above 200 Hz and low-frequency noise below 10 Hz

In [17]: 
# [Signal Filtering]
low_cutoff = 10 # Hz
high_cutoff = 200 # Hz

# Application of the signal to the filter.
signal_filtered = bsnb.bandpass(art_signal, low_cutoff, high_cutoff, order=4, fs = sr)

D.4 - Check the new frequency content of the filtered signal

In [18]: 
#Comparison of Power sprectrum for signal and filtered signal
bsnb.plot_before_after_filter(art_signal, sr, order=4, band_begin=10, band_end=200, x_lim=[0, 500],y_lim=[0, 4e6], show_plot=True,orientation="same");

D.5 - Comparison between the original and filtered signal

In [19]: 
bsnb.plot([list(time), list(time)], [list(art_signal), list(signal_filtered)], legend_label=["Signal", "Filtered"], grid_plot=False, opensignals_style=True, x_axis_label="Time (s)", y_axis_label="RAW Data Samples", x_range=(23, 23.7))

As it can be seen in the power spectrum, most reduction takes place in the frequencies above the cutoff frequency of 200 Hz.

For more information on signal filtering, see Jupyter Notebook entitled "Digital Filtering - A Fundamental Pre-Processing Step" .

E. Event Detection

For the detection of muscular activations, a very simple double-threshold algorithm is defined in the following topics.

E.1 - Signal rectification

In [20]: 
# [Signal Rectification]
rect_signal = absolute(signal_filtered)

E.2 - Envelope Determination

In [21]: 
# Smoothing level [Size of sliding window]
smoothing_level_perc = 20 # Percentage.
smoothing_level = int((smoothing_level_perc / 100) * sr)

#Smooth the signal
signal_smooth = bsnb.smooth(rect_signal, smoothing_level, window='hanning')

E.3 - Double Threshold definition
The first threshold is defined to detect the onset of the muscular activation and the second threshold is defined to detect the offset of the muscular activation of the signal:

Threshold 1 = average + standard deviation

Threshold 2 = average

In [22]: 
# [Threshold]
avg_signal = average(signal_smooth)
std_signal = std(signal_smooth)

thresh_1 = avg_signal + std_signal #Onset of muscular activation
thresh_2 = avg_signal  #Offset of muscular activation
In [23]: 
from bokeh.layouts import layout
from bokeh.plotting import show

fig_list = bsnb.plot([list(time), list(time), list(time)], [list(int((0.5 * max(signal_filtered)) / max(signal_smooth)) * array(signal_smooth)), list(signal_smooth), list(signal_smooth)], title=["Original and Smoothed Signals", "Threshold 1", "Threshold 2"], grid_plot=True, grid_lines=2, grid_columns=2, hor_lines=[[], [thresh_1], [thresh_2]], opensignals_style=True, show_plot=False, x_axis_label="Time (s)", y_axis_label=["Raw Data", "Raw Data", "Raw Data"], get_fig_list=True,x_range=(0,10))
bsnb.opensignals_style(fig_list)

fig_list[0].line(time, signal_filtered)
grid_plot = layout(children=[[fig_list[0]], [fig_list[1], fig_list[2]]], sizing_mode="scale_width")

show(grid_plot)

E.4 - Identification of starting and ending samples of each muscular activation period

In [24]: 
# Initialisation of the list that will store the start and end samples of each muscular activation period.
index_on_off=[]

# Flag that identifies when the onset was detected and the offset searching stage should begin.
flag_onset = False

# Temporary variable that will store the onset until its storage on the definitive list.
temp_onset = None

#returning the index for either onset or offset of activation
for index, sample_value in enumerate(signal_smooth):
    # Identification of the onset index.
    if sample_value > thresh_1 and flag_onset is False:
        # Temporary storage of the onset index.
        temp_onset = index
        
        # Update flag.
        flag_onset = True
    # Identification of the offset index.
    elif sample_value < thresh_2 and flag_onset is True: 
        # Definitive storage of the onset and offset indexes.
        index_on_off.append((temp_onset, index))
        
        # Update flag.
        flag_onset = False

E.5 - Generation of a binary signal intended to identify each muscular activation period

In [25]: 
# Initialisation of the binary signal.
binary_signal = zeros(len(signal_smooth))

# Change the binary signal to 1 for all samples between the onset and offset indexes.
for onset_offset in index_on_off:
    binary_signal[onset_offset[0]:onset_offset[1]] = 1
In [26]: 
bsnb.plot_simple_threshold_algorithm(time, signal_smooth, binary_signal, sr, thresh_1, thresh_2)

The plot illustrates the onset of activation with Threshold 1 and the offset of activation with Threshold 2 .

F. Parameter Extraction

List of EMG analysis parameters:

  • Number of Muscular Activations;
  • Maximum, Minimum and Average duration of muscular activations;
  • Minimum, Maximum, Average and Standard Deviation values of EMG samples;
  • Root Mean Square (RMS) and Area under curve;
  • Total Power, Maximum Frequency and Median Frequency;

F.1 - Detection and accounting of muscular activations

In [27]: 
#Begin and end samples of muscular activation periods
act_begin = array(index_on_off)[:, 0]
act_end = array(index_on_off)[:, 1]

# Time instants where each muscular activation starts and ends.
muscular_activation_begin = array(time)[act_begin]
muscular_activation_end = array(time)[act_end]
In [28]: 
 bsnb.plot([list(time)[::10], list(time)[::10]], [list(array(binary_signal) *max(signal_uv))[::10], list(signal_uv)[::10]], legend=["Activation Signal", "EMG Signal"], grid_plot=False, opensignals_style=True, x_axis_label="Time (s)", y_axis_label="RAW EMG")
In [29]: 
#pip install sty and restart kernel
from sty import fg, rs

# Number of activation periods
print (fg(98,195,238) + "\033[1mNumber of Muscular Activations: \033[0m" + fg.rs + str(len(muscular_activation_begin)))

Number of Muscular Activations: 21

F.2 - Maximum, Minimum and Average duration of muscular activation periods

In [30]: 
# Muscular Activation Duration
muscular_activation_time = muscular_activation_end - muscular_activation_begin

# Parameter extraction
max_time = max(muscular_activation_time)
min_time = min(muscular_activation_time)
avg_time = average(muscular_activation_time)
std_time = std(muscular_activation_time)
In [31]: 
print (fg(98,195,238) + "\033[1m[Maximum, Minimum, Average] duration of Muscular Activations \033[0m" + fg.rs + " = [" + str(max_time) + ", " + str(min_time) + ", " + str(avg_time) + "] s")
print (fg(98,195,238) + "\033[1mStandard Deviation \033[0m" + fg.rs + "= " + str(std_time) + " s")

[Maximum, Minimum, Average] duration of Muscular Activations  = [0.27000308222696745, 0.18600212331190846, 0.21790724943044046] s
Standard Deviation = 0.021921982807397752 s
In [32]: 
bsnb.plot_emg_graphical_durations(max_time, min_time, avg_time, std_time)

F.3 - Maximum, Minimum, Average and Standard Deviation of EMG sample values

In [33]: 
# Maximum
max_sample_value = max(signal_uv)

# Minimum
min_sample_value = min(signal_uv)

# Average and Standard Deviation
avg_sample_value = average(signal_uv)
std_sample_value = std(signal_uv)

time_param_dict = {"Maximum EMG": max_sample_value, "Minimum EMG": min_sample_value, "Average EMG": avg_sample_value, "Standard Deviation EMG": std_sample_value}
In [34]: 
print (fg(98,195,238) + "\033[1m[Maximum, Minimum, Average, Standard Deviation] mV \033[0m" + fg.rs + " = [" + str(max_sample_value) + ", " + str(min_sample_value) + ", " + str(avg_sample_value) + ", " + str(std_sample_value) + "] mV")

[Maximum, Minimum, Average, Standard Deviation] mV  = [1.467041015625, -1.4990386962890625, 0.0019273470525872217, 0.07100350404829539] mV
In [35]: 
bsnb.plot_emg_graphical_statistical(time, signal_uv, max_sample_value, min_sample_value, avg_sample_value, std_sample_value, y_range=[-1.6, 1.6])

F.4 - Root Mean Square and Area under the curve (Signal Intensity Estimators)

In [36]: 
# Root Mean Square
rms = sqrt(sum(signal_uv * signal_uv) / len(signal_uv))

# Area under the curve
area = cumtrapz(signal_uv)
In [37]: 
print (fg(98,195,238) + "\033[1mRoot Mean Square \033[0m" + fg.rs + " = " + str(rms) + " mV")
print (fg(98,195,238) + "\033[1mArea \033[0m" + fg.rs + " = " + str(area[-1]) + " mV.s")

Root Mean Square  = 0.07102965756497363 mV
Area  = 168.83223724365234 mV.s
In [38]: 
from numpy import concatenate

#zoom in plot
bsnb.plot_emg_rms_area(concatenate((time[::10], [time[-1]])), concatenate((signal_uv[::10], [signal_uv[-1]])), rms, area[::10])

F.5 - Total power and some reference points on the frequency domain

In [39]: 
# Signal Power Spectrum
f, P = welch(signal_uv, fs=sr, window='hanning', noverlap=0, nfft=int(256.))

# Total Power and Median Frequency (Frequency that divides the spectrum into two regions with equal power)
area_freq = cumtrapz(P, f, initial=0)
total_power = area_freq[-1]
median_freq = f[where(area_freq >= total_power / 2)[0][0]]
f_max = f[argmax(P)]
In [40]: 
print (fg(98,195,238) + "\033[1mTotal Power \033[0m" + fg.rs + " = " + str(total_power))
print (fg(98,195,238) + "\033[1m[Median Frequency, Maximum Power Frequency] \033[0m" + fg.rs + " = [" + str(median_freq) + ", " + str(f_max) + "] Hz")

Total Power  = 0.004399436828856141
[Median Frequency, Maximum Power Frequency]  = [101.5625, 109.375] Hz
In [41]: 
bsnb.plot_emg_spect_freq(f, P, f_max, median_freq)

This procedure can be automatically done by emg_parameters function in extract module of biosignalsnotebooks package

In [42]: 
bsnb.emg_parameters(signal_uv, sr, raw_to_mv=False)
Out[42]: 
{'Number of Muscular Activations': 21,
 'Maximum Muscular Activation Duration': 0.5040057534903326,
 'Minimum Muscular Activation Duration': 0.3840043836116891,
 'Average Muscular Activation Duration': 0.4651481670534394,
 'Standard Deviation of Muscular Activation Duration': 0.023901429851217293,
 'Maximum Sample Value': 1.467041015625,
 'Minimum Sample Value': -1.4990386962890625,
 'Average Sample Value': 0.0019273470525872217,
 'Standard Deviation Sample Value': 0.07100350404829539,
 'RMS': 0.07102965756497363,
 'Area': 168.83223724365234,
 'Total Power Spect': 0.004399436828856141,
 'Median Frequency': 101.5625,
 'Maximum Power Frequency': 109.375}

This Notebook was a short overview on the processing steps that can be taken into account for EMG signal processing for the extraction of parameters as well as event detection.

For further and more detailed processing methods such as the TKEO operator, please also see the available Jupyter Notebooks entitled "Event Detection - Muscular Activations (EMG)" and for information regarding Fatigue Evaluation of EMG signals "Fatigue Evaluation - Evolution of Median Power Frequency" .

We hope that you have enjoyed this guide. biosignalsnotebooks is an environment in continuous expansion, so don"t stop your journey and learn more with the remaining Notebooks !

**Auxiliary Code Segment (should not be replicated by the user)**

In [43]: 
from biosignalsnotebooks.__notebook_support__ import css_style_apply
css_style_apply()

.................... CSS Style Applied to Jupyter Notebook .........................
Out[43]: